﻿#include "CalculateArea.h"

#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/filter/geometry_refinement/worklet/Triangulate.h>
#include <vtkm/worklet/WorkletMapTopology.h>

struct CalAreaWorklet : public vtkm::worklet::WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn, FieldInPoint, FieldOutCell);
  using ExecutionSignature = void(_1, _2, _3);

  template<typename Fty, typename... Args>
  VTKM_EXEC void CallOnCell(vtkm::CellShapeTagGeneric tag, Fty&& f, Args&&... args) const
  {
    switch (tag.Id)
    {
      case vtkm::CELL_SHAPE_TRIANGLE:
      {
        using CellShapeTag = vtkm::CellShapeTagTriangle;
        //f(CellShapeTag{}, std::forward<Args>(args)...);
      }
      break;
      case vtkm::CELL_SHAPE_QUAD:
      {
        using CellShapeTag = vtkm::CellShapeTagQuad;
        f(CellShapeTag{}, std::forward<Args>(args)...);
      };
      break;
      default:
      {
        //EXPECT_TRUE(false); // 不希望运行到这里
      }
      break;
    }
  }

  template<typename PointVecType, typename CellShapeTag>
  VTKM_EXEC void operator()(CellShapeTag& shape, const PointVecType& pts, Real& area) const
  {
    vtkm::CellShapeTagGeneric shagpe_generic{ shape.Id };
    //std::cout << int(shape.Id) << std::endl;
    //CallOnCell(shagpe_generic, [&area, &pts](auto tag) { area = CellArea(tag, pts); });
    area = CellArea(vtkm::CellShapeTagQuad{}, pts);
  }

  template<typename PointVecType>
  VTKM_EXEC Real CellArea(vtkm::CellShapeTagTriangle&&, const PointVecType& points) const
  {
    auto directed_area = vtkm::Cross(points[0] - points[1], points[2] - points[1]);
    return vtkm::Magnitude(directed_area);
  }

  template<typename PointVecType>
  VTKM_EXEC Real CellArea(vtkm::CellShapeTagQuad&& tag, const PointVecType& points) const
  {
    auto area = vtkm::Magnitude(vtkm::Cross(points[0] - points[1], points[2] - points[1]));
    area += vtkm::Magnitude(vtkm::Cross(points[2] - points[3], points[0] - points[3]));
    return area / 2.;
  }
};


CalculateArea::CalculateArea() {}

void CalculateArea::Init() {}

void CalculateArea::Evolve()
{
  std::vector<int> _dims{ 10000, 10000 };
  std::vector<Real> _x_ranges{ 0, 1 };
  std::vector<Real> _y_ranges{ 0, 1 };
  std::vector<Real> _z_ranges{ 0, 1 };

  Vec2f start{ _x_ranges[0], _y_ranges[0] };
  Vec2f end{ _x_ranges[1], _y_ranges[1] };
  Vec2f space = end - start;
  Id2 point_dim = { _dims[0] + 1, _dims[1] + 1 };

  for (size_t i = 0; i < 2; i++)
    space[i] = (end[i] - start[i]) / _dims[i];

  auto data_set = vtkm::cont::DataSetBuilderUniform::Create(point_dim, start, space);

  vtkm::worklet::Triangulate triangulate;
  auto out_cellset = triangulate.Run(data_set.GetCellSet());

  ArrayHandle<Real> area;
  std::cout << "单元数：" << data_set.GetCellSet().GetNumberOfCells() << std::endl;
  std::cout << "节点数：" << data_set.GetCoordinateSystem().GetNumberOfValues() << std::endl;

  vtkm::cont::Timer timer;
  timer.Start();
  vtkm::cont::Invoker{}(CalAreaWorklet{},
                        data_set.GetCellSet(),
                        data_set.GetCoordinateSystem().GetDataAsMultiplexer(),
                        area);

  auto toal_area = vtkm::cont::Algorithm::Reduce(area, 0.0);

  std::cout << "时间开销：" << timer.GetElapsedTime() << std::endl;
  std::cout << "total area：" << toal_area << std::endl;
}
